Illustration of iterative linear solver behavior on simple 1D and 2D problems

In geometry processing, numerical optimization methods often involve solving sparse linear systems of equations. These linear systems have a structure that strongly resembles to adjacency graphs of the underlying mesh. We observe how classic linear solvers behave on this specific type of problems. For the sake of simplicity, we minimise either the squared gradient or the squared Laplacian, evaluated by finite differences on a regular 1D or 2D grid. We observed the evolution of the solution for both energies, in 1D and 2D, and with different solvers: Jacobi, Gauss-Seidel, SSOR (Symmetric successive over-relaxation) and CG (conjugate gradient [She94]). Plotting results at different iterations allows to have an intuition of the behavior of these classic solvers. Introduction When facing an optimization problem, the simplest approach is often to iterate, locally improving current solution. For example, heat diffusion can be done by iteratively replacing each value by the average value of its neighbors. With more background on numerical optimization, one would like to formulate this optimization as a linear problem. With such abstraction, it is possible to use solvers as black boxes: we do not have an interpretation of the evolution of the solution during the iterations. We are interested in observing the behavior of basic linear solvers for very simple cases. It does not necessarily impact the way to use them, but provides some intuition about the magic of linear solvers in the context of minimizing energies on a mesh. Our meshes are regular 1D and 2D grids of size n, and variables are attached to vertices. We address the coordinates of the unknown vector x as follows: in 1D xi is the i th element of x and in 2D, the value xi,j of the vertex located at the i th line and k column is the (n ∗ i + j) element of x. Our tests are performed on the minimization of two energies: • Gradient energy: The sum of squared difference between adjacent samples: Σi(xi − xi+1) in 1D and Σi,j(xi,j − xi+1,j) + Σi,j(xi,j − xi,j+1) in 2D. • Laplacian energy: The sum of squared difference between a sample and the average of its neighbors: Σi(xi−xi−1/2−xi+1/2) in 1D. In 2D, it is the same on energy for boundary vertices, and Σi,j(xi,j − xi−1,j/2 − xi+1,j/2) + (xi,j − xi,j−1/2 − xi,j+1/2) over other vertices. Note that for the Laplacian energy with the right number of constraints (2 in 1D, and 3 in 2D), the minimization of the first energy is equivalent to the minimization of the second energy. 1 1D problems Observing 1D problems allows to visualize the convergence by a 3D surface. In all figures in this section, the first axis (horizontal) represents the domain — a 1D grid, the second axis (vertical) represents the value associated to each vertex, and the third axis (depth) is the time. 1 ar X iv :1 51 0. 01 11 8v 1 [ cs .N A ] 5 O ct 2 01 5 1.1 Gradient energy The 1D grid has n vertices, here we develop the matrices for n = 5. The problem is formalized by this set of equations Ax = b :  0 −1 0 0 0 0 1 −1 0 0 0 0 1 −1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 x =  −2 0 0 6 2 6  The 4 first lines represent the objective. To enforce boundary constraints x1 = 2 and x5 = 6, we replaced coefficients applied to x1 and x5 by a contribution to the right hand side. We also add the last two lines to explicitly set x1 = 2 and x5 = 6. In the least squares sense, it produces the new system to solve A >Ax = A>b i.e.  1 0 0 0 0 0 2 −1 0 0 0 −1 2 −1 0 0 0 −1 2 0 0 0 0 0 1 x =  2 2 0 6 6  The condition number of A>A is 5.9 for n = 5, and 1000 for n = 50. These numbers are not high enough to produce numerical instabilities. This first experience (Fig. 1) shows the behavior of each solver with the same matrix. The boundary conditions (locked coordinates of x) are set in the way that the solution is the straight line defined by xi = i(xn − x1)/n. We observe the expected relative speed of convergence i.e. Jacobi < Gauss-Seidel < SSOR < Conjugate Gradient. We also visualize the impact of the order of the coordinates in Gauss-Seidel and SSOR: when only the left side is constrained (locked x1), the first iteration propagates the constraints to all coordinates whereas constraining only xn take n iteration to affect all the variables. 1.2 Laplacian energy As for the gradient energy case, the 1D grid has n vertices, and we develop the matrices for n = 5. The problem is formalized by this set of equations Ax = b :  0 2 −1 0 0 0 −1 2 −1 0 0 0 −1 2 0 1 0 0 0 0 0 0 0 0 1 x =  −2 0 6 2 6  The 3 first lines represent the objective. As for the “gradient energy”, we replaced coefficients applied to x1 and x5 by a contribution to the right hand side to enforce x1 = 2 and x5 = 6. We also add the last two lines to explicitly set x1 = 2 and x5 = 6. In least squares sense, it produces the new system to solve A>Ax = A>b i.e.  1 0 0 0 0 0 5 −4 1 0 0 −4 6 −4 0 0 1 −4 5 0 0 0 0 0 1 x =  2 4 −8 12 6  The condition number of A>A is 34.4 for n = 5, and 920000 for n = 50 coordinates. Such numbers produce numerical instabilities slowing down the solvers. In Fig. 3, we observe the same expected relative speed of convergence i.e. Jacobi < Gauss-Seidel < SSOR < Conjugate Gradient. Oscillations of the solution with Gauss-Seidel and SSOR are interesting to observe: it converges faster with higher local oscillations (in the solution space, not in the time dimension). Another interesting observation is that CG required more than 20 iterations to converge. This was unexpected because, while CG is mostly used as an iterative solver, it is also a direct solver that is supposed